SUBROUTINE read_ini_data
  USE Commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'

  REAL(8), ALLOCATABLE :: nonpardisnoctemp(:,:), nonpardistemaxnorcnosav(:,:), nonpardistemaxnorcnosavtemp(:,:), temp(:), initdata_temp(:,:), work(:)
  REAL(8)              :: ptemp1, ptemp2, ptemp3, datatemp(Nmomdata,1), factor, temp1, temp2(Nmomdata), prcont(Nrcont)
  INTEGER              :: i, j, iwealthc, ipr, ifu, irc, lwork, info
  CHARACTER          :: uplo='U'

  !---------------------------------------------------------------
  !   INITIAL CONDITIONS
  !   read from a text file with the following format:
  !1. Municipal id
  !2. Audit (1/0)
  !3. Mayor's age (1 unit corresponds to a term, which is 4 years)
  !4. Mayor's educ (1,2,3)
  !5. Private GDP (industry index) 97-00
  !6. Public funds 97-00
  !7. Mayor's wage
  !8. Number of terms
  !9. Outcome election
  !10. Population thresholds (1,2,3,4,5,6)
  !11. Population dummy for 0-10,000
  !12. Population dummy for 10,000-50,000
  !13. Population dummy for 50,000-100,000
  !14. Public funds 01-04
  !15. Mayor's age 01-04(1 unit corresponds to a term, which is 4 years)
  !16. Mayor's educ 01-04(1,2,3)
  !17. Private GDP (industry index) 01-04
  !18. Candidate's wealth
  !Not included 14. Population dummy for 100,000-300,000
  !Not included 15. Population dummy for 300,000-500,000
  !Not included 16. Population dummy for >500,000
  !---------------------------------------------------------------

  ALLOCATE(countNnterms_mun1(Nnterms), countNnterms_mun2(Nnterms), countNnterms_mun3(Nnterms))

  ALLOCATE(initdata(Nactmun,Nvarinit),initdata_temp(Nactmun,Nvarinit),temp(Nvarinit))
  initdata = 0.0d0

!!$  OPEN(10, file='/u/home/m/mmazzoc/project/MyPapers/Mayors/Estimation/Stata/initial_new_last.raw')
  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/initial_Feb_2020.raw')
  DO i = 1, Nactmun
     READ(10, *) temp(1:Nvarinit)
     initdata(i,1:Nvarinit) = temp(1:Nvarinit)
!!$     IF (myrank==0) WRITE(*,*) 'initdata, i', initdata(i,1:Nvarinit), i
  END DO
  DEALLOCATE(temp)
  CLOSE(10)

  IF (flag_bootstrap == 0) THEN
     ! Average variables to be used when we control for observatble heterogeneity
     av_saving = (SUM(initdata(:,18))/DBLE(Nactmun))/norm
     IF (myrank==0) WRITE(*,'(a,1f20.10)') 'av_saving', av_saving
     av_zzpr = (SUM(initdata(:,5))/DBLE(Nactmun))
     IF (myrank==0) WRITE(*,'(a,1f20.10)') 'av_zzpr', av_zzpr
     av_funds = (SUM(initdata(:,6))/DBLE(Nactmun))/norm
     IF (myrank==0) WRITE(*,'(a,1f20.10)') 'av_funds', av_funds     
  END IF
  
  !---------------------------------------------------------------
  !   DATA USED TO CONSTRUCT THE MATRIX THAT REPRESENTS THE NONPARAMETERIC DISTRIBUTION FOR AGE, EDUC, ETC.
  !   read from a text file with the following format:
  !1. Municipal id
  !2. Mayor's age (1 unit corresponds to a term, which is 4 years)
  !3. Mayor's educ (1,2,3)
  !4. Public funds
  !5. Private GDP
  !6. Number of terms the mayor was in power
  !7. Mayor's wealth
  !---------------------------------------------------------------

  ! SMALL MUNICIPALITIES: Open the data set for the non-parameteric distribution

  ALLOCATE(nonpardistnoc_mun1(Ntotmunnoc_mun1,Nnterms,Nvarnonparnoc))
  nonpardistnoc_mun1 = 0d0

  ALLOCATE(nonpardistnoc_uncond_mun1(Ntotmunnoc_mun1,2))
  nonpardistnoc_uncond_mun1 = 0d0

  ALLOCATE(nonpardisnoctemp(Ntotmunnoc_mun1,Nvarnonparnoctemp),temp(Nvarnonparnoctemp))
  nonpardisnoctemp = 0d0

  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/probdata_1_2019.raw')
  DO i = 1, Ntotmunnoc_mun1
     READ(10, *) temp(1:Nvarnonparnoctemp)
     nonpardisnoctemp(i,1:Nvarnonparnoctemp) = temp(1:Nvarnonparnoctemp)
!!$     IF (myrank==0) WRITE(*,'(a,7f20.10,1i4)') 'nonpardisnoctemp, i', nonpardisnoctemp(i,1:Nvarnonparnoctemp), i
  END DO
  CLOSE(10)

  countNnterms_mun1 = 0
  DO i = 1, Ntotmunnoc_mun1
     DO j = 1, Nnterms
        IF (INT(nonpardisnoctemp(i,6)) == j) THEN
           countNnterms_mun1(j) = countNnterms_mun1(j) + 1
           nonpardistnoc_mun1(countNnterms_mun1(j),j,1) = nonpardisnoctemp(i,2)
           nonpardistnoc_mun1(countNnterms_mun1(j),j,2) = nonpardisnoctemp(i,3)
           nonpardistnoc_mun1(countNnterms_mun1(j),j,3) = nonpardisnoctemp(i,4)/norm
           nonpardistnoc_mun1(countNnterms_mun1(j),j,4) = nonpardisnoctemp(i,5)
           nonpardistnoc_mun1(countNnterms_mun1(j),j,5) = nonpardisnoctemp(i,7)/norm
!!$           IF (myrank==0) WRITE(*,*) 'nonpardist(i,j,:),i,j', countNnterms(j), j, nonpardist(countNnterms(j),j,:)
        END IF
     END DO
  END DO

  !Funds and private GDP are municipality characteristics and not mayor's characteristics; we should therefore not condition on nterms
  nonpardistnoc_uncond_mun1(:,1) = nonpardisnoctemp(:,4)/norm
  nonpardistnoc_uncond_mun1(:,2) = nonpardisnoctemp(:,5)

  DEALLOCATE(nonpardisnoctemp,temp)

  ! MEDIUM MUNICIPALITIES: Open the data set for the non-parameteric distribution

  ALLOCATE(nonpardistnoc_mun2(Ntotmunnoc_mun2,Nnterms,Nvarnonparnoc))
  nonpardistnoc_mun2 = 0d0

  ALLOCATE(nonpardistnoc_uncond_mun2(Ntotmunnoc_mun2,2))
  nonpardistnoc_uncond_mun2 = 0d0

  ALLOCATE(nonpardisnoctemp(Ntotmunnoc_mun2,Nvarnonparnoctemp),temp(Nvarnonparnoctemp))
  nonpardisnoctemp = 0d0

  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/probdata_2_2019.raw')
  DO i = 1, Ntotmunnoc_mun2
     READ(10, *) temp(1:Nvarnonparnoctemp)
     nonpardisnoctemp(i,1:Nvarnonparnoctemp) = temp(1:Nvarnonparnoctemp)
!!$     IF (myrank==0) WRITE(*,*) 'nonpardisnoctemp, i', nonpardisnoctemp(i,1:Nvarnonparnoc), i
  END DO
  CLOSE(10)

  countNnterms_mun2 = 0
  DO i = 1, Ntotmunnoc_mun2
     DO j = 1, Nnterms
        IF (INT(nonpardisnoctemp(i,6)) == j) THEN
           countNnterms_mun2(j) = countNnterms_mun2(j) + 1
           nonpardistnoc_mun2(countNnterms_mun2(j),j,1) = nonpardisnoctemp(i,2)
           nonpardistnoc_mun2(countNnterms_mun2(j),j,2) = nonpardisnoctemp(i,3)
           nonpardistnoc_mun2(countNnterms_mun2(j),j,3) = nonpardisnoctemp(i,4)/norm
           nonpardistnoc_mun2(countNnterms_mun2(j),j,4) = nonpardisnoctemp(i,5)
           nonpardistnoc_mun2(countNnterms_mun2(j),j,5) = nonpardisnoctemp(i,7)/norm
!!$           IF (myrank==0) WRITE(*,*) 'nonpardist(i,j,:),i,j', countNnterms(j), j, nonpardist(countNnterms(j),j,:)
        END IF
     END DO
  END DO

  !Funds and private GDP are municipality characteristics and not mayor's characteristics; we should therefore not condition on nterms
  nonpardistnoc_uncond_mun2(:,1) = nonpardisnoctemp(:,4)/norm
  nonpardistnoc_uncond_mun2(:,2) = nonpardisnoctemp(:,5)

  DEALLOCATE(nonpardisnoctemp,temp)

  ! LARGE MUNICIPALITIES: Open the data set for the non-parameteric distribution

  ALLOCATE(nonpardistnoc_mun3(Ntotmunnoc_mun3,Nnterms,Nvarnonparnoc))
  nonpardistnoc_mun3 = 0d0

  ALLOCATE(nonpardistnoc_uncond_mun3(Ntotmunnoc_mun3,2))
  nonpardistnoc_uncond_mun3 = 0d0

  ALLOCATE(nonpardisnoctemp(Ntotmunnoc_mun3,Nvarnonparnoctemp),temp(Nvarnonparnoctemp))
  nonpardisnoctemp = 0d0

  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/probdata_3_2019.raw')
  DO i = 1, Ntotmunnoc_mun3
     READ(10, *) temp(1:Nvarnonparnoctemp)
     nonpardisnoctemp(i,1:Nvarnonparnoctemp) = temp(1:Nvarnonparnoctemp)
!!$     IF (myrank==0) WRITE(*,*) 'nonpardisnoctemp, i', nonpardisnoctemp(i,1:Nvarnonparnoc), i
  END DO
  CLOSE(10)

  countNnterms_mun3 = 0
  DO i = 1, Ntotmunnoc_mun3
     DO j = 1, Nnterms
        IF (INT(nonpardisnoctemp(i,6)) == j) THEN
           countNnterms_mun3(j) = countNnterms_mun3(j) + 1
           nonpardistnoc_mun3(countNnterms_mun3(j),j,1) = nonpardisnoctemp(i,2)
           nonpardistnoc_mun3(countNnterms_mun3(j),j,2) = nonpardisnoctemp(i,3)
           nonpardistnoc_mun3(countNnterms_mun3(j),j,3) = nonpardisnoctemp(i,4)/norm
           nonpardistnoc_mun3(countNnterms_mun3(j),j,4) = nonpardisnoctemp(i,5)
           nonpardistnoc_mun3(countNnterms_mun3(j),j,5) = nonpardisnoctemp(i,7)/norm
!!$           IF (myrank==0) WRITE(*,*) 'nonpardist(i,j,:),i,j', countNnterms(j), j, nonpardist(countNnterms(j),j,:)
        END IF
     END DO
  END DO

  !Funds and private GDP are municipality characteristics and not mayor's characteristics; we should therefore not condition on nterms
  nonpardistnoc_uncond_mun3(:,1) = nonpardisnoctemp(:,4)/norm
  nonpardistnoc_uncond_mun3(:,2) = nonpardisnoctemp(:,5)

  DEALLOCATE(nonpardisnoctemp,temp)

  !----------------------------------------------------------------------------------
  !   READING IN THE NONPARAMETRIC PROBABILITIES OF CHALLENGER AGE AND EDUCATION
  !   read from a text file with the following format:(specific file for Retage=8)
  !1. probability mayor age = x AND educ = 1 where x=1,2,3
  !   probability mayor age = x AND educ = 2 where x=1,2,3
  !---------------------------------------------------------------
  ALLOCATE(nonpardistageduc(Nage,Neduc))
  nonpardistageduc=0.0d0
  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/candageduc.raw')
  DO i=1,Neduc
     DO j = 1,Nage
        READ(10,*) temp1
        nonpardistageduc(j,i) = temp1
     END DO
  END DO
  CLOSE(10)

  !---------------------------------------------------------------
  !   DATA USED TO CONSTRUCT THE MATRIX THAT REPRESENTS THE NONPARAMETERIC DISTRIBUTION FOR RELATIVE CONTR.
  !   read from a text file with the following format:
  !1. Relative contributions
  !---------------------------------------------------------------

  ALLOCATE(nonpardistc(Ntotmunc))
  nonpardistc = 0.0d0

  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/probcontrib_new_2019.raw')
  DO i = 1, Ntotmunc
     READ(10, *) temp1
     nonpardistc(i) = temp1
!!$     IF (myrank==0) WRITE(*,*) 'nonpardisc, i', nonpardisc(i), i
  END DO
  CLOSE(10)

!!$     !Count how many observations we have with Nnterms = i
!!$     DO i = 1, Nnterms
!!$        if (myrank == 0) write(*,*) 'countNnterms(i), 1, 2, 3', i, countNnterms_mun1(i), countNnterms_mun2(i), countNnterms_mun3(i)
!!$     END DO

  !---------------------------------------------------------------
  !   NONPARAMETERIC DISTRIBUTION FOR RELATIVE CONTRIBUTIONS, FUNDS, AND PRIVATE GDP
  !   read from a text file with the following format:
  !1. Relative Contributions (1-3)
  !2. Private GDP (1-3)
  !3. Public funds (1-3)
  !4. Probabilities
  !---------------------------------------------------------------

  !WE HAVE TO COMPUTE THE CORRECT PROBABILITY FOR pwealth
  pwealth = 1d0/DBLE(Nwealth)
  !The following is the correct prcont since we use quartiles for the grid
  prcont = 1d0/DBLE(Nrcont)

  ALLOCATE(nonpardistemax(Nwealth,Nfunds,Npr,Nrcont), nonpardistemaxnosav(Nfunds,Npr,Nrcont), &
           nonpardistemaxnorcnosav(Nfunds,Npr), nonpardistemaxnorc(Nwealth,Nfunds,Npr))
  nonpardistemax = 0d0
  nonpardistemaxnorcnosav = 0d0
  nonpardistemaxnorc = 0d0

  !READ IN THE NONPARAMETRIC DISTRIBUTION COMPUTED USING THE DATA
  OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/probabilities_new_last.raw')
  ptemp1 = 0d0
  ptemp2 = 0d0
  ptemp3 = 0d0

!!$  DO irc = 1, Nrcont
!!$     DO ifu = 1, Nfunds
!!$        DO ipr = 1, Npr
!!$           READ(10, *) nonpardistemaxnosav(ifu,ipr,irc)
!!$           ptemp1 = ptemp1 + nonpardistemaxnosav(ifu,ipr,irc)
!!$
!!$           DO iwealthc = 1, Nwealth
!!$              nonpardistemax(iwealthc,ifu,ipr,irc) = nonpardistemaxnosav(ifu,ipr,irc)*pwealth(iwealthc)
!!$
!!$              ptemp2 = ptemp2 + nonpardistemax(iwealthc,ifu,ipr,irc)
!!$
!!$              !IF (myrank==0) WRITE(*,*) 'nonpardistemaxnosav,nonpardistemax',nonpardistemaxnosav(ifu,ipr,irc),nonpardistemax(iwealthc,ifu,ipr,irc)
!!$
!!$           END DO
!!$        END DO
!!$     END DO
!!$  END DO
!!$  CLOSE(10)
!!$
!!$  IF (ptemp1 < 0.99d0 .OR. ptemp1 > 1.01d0 .OR. ptemp2 < 0.99d0 .OR. ptemp2 > 1.01d0) THEN
!!$     IF (myrank == 0) WRITE(*,*) 'pemaxnosav, pemax', ptemp1, ptemp2
!!$  END IF

  IF (Npr > 1) THEN 
     DO ifu = 1, Nfunds
        DO ipr = 1, Npr
           READ(10,*) temp1
           nonpardistemaxnorcnosav(ifu,ipr) = temp1
        END DO
     END DO
  ELSE
     ALLOCATE(nonpardistemaxnorcnosavtemp(Nfunds,3))
     nonpardistemaxnorcnosavtemp(Nfunds,3) = 0d0
     DO ifu = 1, Nfunds
        ! In the data, we use Npr = 3
        DO ipr = 1, 3
           READ(10,*) temp1
           nonpardistemaxnorcnosavtemp(ifu,ipr) = temp1
        END DO
        nonpardistemaxnorcnosav(ifu,1) = SUM(nonpardistemaxnorcnosavtemp(ifu,:))
     END DO
  END IF

  DO ifu = 1, Nfunds
     DO ipr = 1, Npr

        ptemp1 = ptemp1 + nonpardistemaxnorcnosav(ifu,ipr)

        DO irc = 1, Nrcont

           nonpardistemaxnosav(ifu,ipr,irc) = nonpardistemaxnorcnosav(ifu,ipr)*prcont(irc)

           ptemp2 = ptemp2 + nonpardistemaxnosav(ifu,ipr,irc)

           DO iwealthc = 1, Nwealth

              nonpardistemax(iwealthc,ifu,ipr,irc) = nonpardistemaxnosav(ifu,ipr,irc)*pwealth(iwealthc)

              nonpardistemaxnorc(iwealthc,ifu,ipr) = nonpardistemaxnorcnosav(ifu,ipr)*pwealth(iwealthc)

              ptemp3 = ptemp3 + nonpardistemax(iwealthc,ifu,ipr,irc)

              !IF (myrank==0) WRITE(*,*) 'nonpardistemaxnosav,nonpardistemax',nonpardistemaxnosav(ifu,ipr,irc),nonpardistemax(iwealthc,ifu,ipr,irc)

           END DO
        END DO
     END DO
  END DO
  CLOSE(10)

!!$  IF (ptemp1 < 0.99d0 .OR. ptemp1 > 1.01d0 .OR. ptemp2 < 0.99d0 .OR. ptemp2 > 1.01d0 .OR. ptemp3 < 0.99d0 .OR. ptemp3 > 1.01d0) THEN
!!$     IF (myrank == 0) WRITE(*,*) 'pemax, pemaxnosav, pemaxnorcnosav', ptemp1, ptemp2, ptemp3
!!$  END IF

  ALLOCATE(nonpardist_fu_pr(Nfunds,Npr),nonpardist_wealth_fu_pr(Nwealth,Nfunds,npr))
  !We construct some of the arrays used in probemax.f90 and probemaxtemp.f90
  DO ipr = 1, Npr
     DO ifu = 1, Nfunds
        nonpardist_fu_pr(ifu,ipr) = SUM(nonpardistemax(:,ifu,ipr,:))
!!$        nonpardist_fu_pr(ifu,ipr) = SUM(nonpardistemaxnosav(ifu,ipr,:))
     END DO
  END DO

  DO ipr = 1, Npr
     DO ifu = 1, Nfunds
        DO iwealthc = 1, Nwealth
           nonpardist_wealth_fu_pr(iwealthc,ifu,ipr) = SUM(nonpardistemax(iwealthc,ifu,ipr,:))
        END DO
     END DO
  END DO

  !If using probit for electoral rule
!!$  !Read the moment conditions in the data
!!$  OPEN(20, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/moments_lastt.raw')
!!$  DO j=1, Nmomdata
!!$     READ(20,*) data(j,1)
!!$  END DO
!!$  CLOSE(20) 
!!$  
!!$  !Read covariance matrix of the moment conditions in the data
!!$  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_lastt.raw')
!!$  DO j=1, Nmomdata
!!$     READ(25,*) dataV(j,:)
!!$  END DO
!!$  CLOSE(25) 

  !If using a linear specification for the electoral rule
  !Read the moment conditions in the data
!!$  OPEN(20, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/moments_linear.raw')
  OPEN(20, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/moments_Aug_2023.raw')
  DO j=1, Nmomdata
     READ(20,*) temp1
     data(j,1) = temp1
  END DO
  CLOSE(20) 

  IF (flag_bootstrap == 1) THEN

     OPEN(10, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/moments_bootstrap.raw')
     data(:,1) = 0d0
     DO i = 1, nboot
        READ(10,*) data(:,1)
        IF (myrank==0) write(*,'(a,1i5,11f20.10)') 'moments, bootstrap 0', nboot, data(:,1)
     END DO
     CLOSE(10)

  END IF

  
  !Read covariance matrix of the moment conditions in the data
!!$  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_linear.raw')
!!$  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_probit_June2019.raw')
  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_Aug_2023.raw')
  DO j=1, Nmomdata
     READ(25,*) temp2(:)
     dataV(j,:) = temp2(:)
!!$     IF (myrank == 0) WRITE(*,'(a,12f20.10)') 'dataV', dataV(j,:)
  END DO
  CLOSE(25) 

  datatemp(:,1) = data(:,1)
  DO j=1, Nmomdata
     factor = 1d0
     DO WHILE (ABS(datatemp(j,1)) > 1d0)
        factor = factor * 0.1d0
        datatemp(j,1) = data(j,1)*factor
!!$        IF (myrank == 0) WRITE(*,'(a,2f20.10)') 'factor', factor, datatemp(j,1)
     END DO
     norm_factor(j) = factor
       IF (flag_bootstrap == 0 .AND. myrank == 0) WRITE(*,'(a,2f20.10)') 'norm_factor', norm_factor(j), data(j,1)
  END DO

  !stealing 2nd term
  norm_factor(1) = norm_factor(1)*10d0
  norm_factor(2) = norm_factor(2)*10d0
!!$  norm_factor(3) = norm_factor(3)*10d0
!!$  norm_factor(7) = norm_factor(7)*10d0
  norm_factor(8) = norm_factor(8)*10d0
!!$  norm_factor(10) = norm_factor(10)*100d0

  ! We compute the inverse of the variance-covariance matrix to use as efficient weighting matrix
  V = dataV(1:Nmom,1:Nmom)
  INV_V = 0d0
  DO j=1, Nmom
     INV_V(j,j)=1d0
  END DO

  IF (flag_eff_weight == 1) THEN
     !IF WE USE THE EFFICIENT MATRIX
     LWORK = Nmom
     IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
     ALLOCATE(WORK(LWORK))
     CALL DPOSV(uplo,Nmom,Nmom,V,Nmom,INV_V,Nmom,info)
     IF (info /= 0) THEN
        write(*,*) 'WEIGHT MATRIX COULD NOT INVERT, info', info
        IF(istop == 1) STOP
     END IF
  END IF

END SUBROUTINE read_ini_data
